###################################################
## SonPark_replication.R
## Replication R code for manuscript
## Date: "Thu Aug 11 10:45:26 2022"
###################################################
rm(list=ls())
start_time <- Sys.time()
###################################################
## Load packages
###################################################
library(ggplot2)
library(tibble)
library(gridExtra)
library(readxl)
library(dplyr)
library(RColorBrewer)
library(ggthemes)
library(kableExtra)
library(broom)
library(jtools)
library(tidyr)
library(car)
library(ggthemes)
library(stargazer)
library(xtable)
library(ggh4x)
library(gridExtra)
library(visreg)
library(dplyr)
library(arsenal) 
library(table1)
library(Rmisc)
library(sjPlot)

###################################################
## Load data
###################################################
load("SonPark_Replication.RData")

## Treatment level information
Treatment.level.2019 <- c("Conditional Military Punishment",
                     "Reduced Nuclear Threat",
                     "Economic Sanction",
                     "Resolution on ROK's NPT Violation",
                     "Nuclear Technology Sanction",
                     "Elite Opposition to Proliferation",
                     "Public Opposition to Proliferation")
df2019$Treatment <- factor(df2019$Treatment, levels=Treatment.level.2019)
df2019$Treatment <- relevel(df2019$Treatment, ref = "Public Opposition to Proliferation")

Treatment.level.2022 = c("Conditional Military Punishment",
                           "Reduced Nuclear Threat",
                           "Economic Sanction",
                           "Resolution on ROK's NPT Violation",
                           "Nuclear Technology Sanction",
                           "Elite Opposition to Proliferation",
                           "Public Opposition to Proliferation",
                           "Control")

df2022$Treatment <- factor(df2022$Treatment, levels=Treatment.level.2022)
df2022$Treatment <- relevel(df2022$Treatment, ref = "Control")


###################################################
## Table 3
###################################################
df.tab3 <- df2019%>% group_by(Treatment) %>%
    dplyr::summarise(N = n(),
                     opinion.change.mean = mean(opinion.change),
                     opinion.change.se = sd(opinion.change)/sqrt(N),
                     behavior.change.mean = mean(behavior.change),
                     behavior.change.se = sd(behavior.change)/sqrt(N)
                     )
df.tab3.out <- t(df.tab3[, -1])[-1,]
df.tab3.out <- rbind(df.tab3.out, df.tab3$N)
colnames(df.tab3.out) <- df.tab3$Treatment
xtable(df.tab3.out, digits=2)

###################################################
## Table 4
###################################################
#### (1) logistic regression on opinion change ####
fit1 <- glm(opinion.change ~ Treatment, data=df2019, family = "binomial") ;summary(fit1)
fit2 <- glm(opinion.change ~ Treatment + age + male + ideology  + political.knowledge +
              income + education, data=df2019, family = "binomial") ;summary(fit2)

#### (2) logistic regression on behavior change ####
fit3 <- glm(behavior.change ~ Treatment, data=df2019, family = "binomial") ;summary(fit3)
fit4 <- glm(behavior.change ~ Treatment+ age + male + ideology  + political.knowledge +
              income + education, data=df2019, family = "binomial") ;summary(fit4)
fit5 <- glm(behavior.change ~ Treatment + relevel(choice.at.treat, ref="No Change") + age+ male + ideology  + political.knowledge +
              income + education, data=df2019, family = "binomial") ;summary(fit5)

#### (3) logistic regression on attitude change ####
fit6 <- glm(attitude.change ~ Treatment, data=df2019, family = "binomial") ;summary(fit6)
fit7 <- glm(attitude.change ~ Treatment + age + male + ideology  + political.knowledge +
              income + education, data=df2019, family = "binomial") ;summary(fit7)

stargazer(fit1, fit2, fit3, fit4, fit5, fit6, fit7)


###################################################
## Table 5
###################################################

tab.mat <- matrix(NA, 3, 2)
######## 2019 Experiemnt ######## 
tab.mat[1,1] <- round(table(df2019$opinion.change)['1']/sum(table(df2019$opinion.change))*100, 2)
tab.mat[2,1] <- round(table(df2019$behavior.change)['1']/sum(table(df2019$behavior.change))*100, 2)
tab.mat[3,1] <- round(table(df2019$attitude.change)['1']/sum(table(df2019$attitude.change))*100, 2)

######## 2022 Experiemnt ########
df2022.nocontrol <- df2022 %>% filter(Treatment!="Control") ## Data without Control group 
tab.mat[1,2] <- round(table(df2022.nocontrol$opinion.change)['1']/sum(table(df2022.nocontrol$opinion.change))*100, 2)
tab.mat[2,2] <- round(table(df2022.nocontrol$behavior.change)['1']/sum(table(df2022.nocontrol$behavior.change))*100, 2)
tab.mat[3,2] <- round(table(df2022.nocontrol$attitude.change)['1']/sum(table(df2022.nocontrol$attitude.change))*100, 2)

xtable(tab.mat)

###################################################
## Table 6
###################################################
Changer <- df2019 %>% filter(df2019$opinion.change=='1')
Non_Changer <- df2019 %>% filter(df2019$opinion.change=='0')


tab6.mat <- matrix(NA, 14, 3)

tab6.mat[1, 1] <- mean(Changer$DRPK.Threat) 
tab6.mat[1, 2] <- mean(Non_Changer$DRPK.Threat) 
tp <- t.test(Changer$DRPK.Threat,Non_Changer$DRPK.Threat) 
tab6.mat[1, 3] <- ifelse(tp$p.value < 0.05, "Yes", "No")


tab6.mat[2, 1] <- mean(Changer$US.Credibility) 
tab6.mat[2, 2] <- mean(Non_Changer$US.Credibility) 
tp <- t.test(Changer$US.Credibility,Non_Changer$US.Credibility) 
tab6.mat[2, 3] <- ifelse(tp$p.value < 0.05, "Yes", "No")

tab6.mat[3, 1] <- mean(Changer$national.security) ## 4.289428
tab6.mat[3, 2] <- mean(Non_Changer$national.security) ##4.584906
tp <- t.test(Changer$national.security,Non_Changer$national.security) ##p-value = 9.537e-14
tab6.mat[3, 3] <- ifelse(tp$p.value < 0.05, "Yes", "No")

tab6.mat[4, 1] <- mean(Changer$nuclear.deterrence) 
tab6.mat[4, 2] <- mean(Non_Changer$nuclear.deterrence) 
tp <- t.test(Changer$nuclear.deterrence,Non_Changer$nuclear.deterrence) 
tab6.mat[4, 3] <- ifelse(tp$p.value < 0.05, "Yes", "No")

tab6.mat[5, 1] <- mean(Changer$self.defense) 
tab6.mat[5, 2] <- mean(Non_Changer$self.defense) 
tp <- t.test(Changer$self.defense,Non_Changer$self.defense) 
tab6.mat[5, 3] <- ifelse(tp$p.value < 0.05, "Yes", "No")

tab6.mat[6, 1] <- mean(Changer$diplomatic.leverage) 
tab6.mat[6, 2] <- mean(Non_Changer$diplomatic.leverage) 
tp <- t.test(Changer$diplomatic.leverage,Non_Changer$diplomatic.leverage) 
tab6.mat[6, 3] <- ifelse(tp$p.value < 0.05, "Yes", "No")

tab6.mat[7, 1] <- mean(Changer$national.pride) 
tab6.mat[7, 2] <- mean(Non_Changer$national.pride) 
tp <- t.test(Changer$national.pride,Non_Changer$national.pride) 
tab6.mat[7, 3] <- ifelse(tp$p.value < 0.05, "Yes", "No")

tab6.mat[8, 1] <- mean(Changer$international.status) 
tab6.mat[8, 2] <- mean(Non_Changer$international.status) 
tp <- t.test(Changer$international.status,Non_Changer$international.status) 
tab6.mat[8, 3] <- ifelse(tp$p.value < 0.05, "Yes", "No")

tab6.mat[9, 1] <- mean(Changer$NPT.Legality) 
tab6.mat[9, 2] <- mean(Non_Changer$NPT.Legality) 
tp <- t.test(Changer$NPT.Legality,Non_Changer$NPT.Legality) 
tab6.mat[9, 3] <- ifelse(tp$p.value < 0.05, "Yes", "No")

tab6.mat[10, 1] <- mean(Changer$NPT.Sanction) 
tab6.mat[10, 2] <- mean(Non_Changer$NPT.Sanction) 
tp <- t.test(Changer$NPT.Sanction,Non_Changer$NPT.Sanction) 
tab6.mat[10, 3] <- ifelse(tp$p.value < 0.05, "Yes", "No")

tab6.mat[11, 1] <- mean(Changer$NPT.Reputation) 
tab6.mat[11, 2] <- mean(Non_Changer$NPT.Reputation) 
tp <- t.test(Changer$NPT.Reputation,Non_Changer$NPT.Reputation) 
tab6.mat[11, 3] <- ifelse(tp$p.value < 0.05, "Yes", "No")

tab6.mat[12, 1] <- mean(Changer$Morality) 
tab6.mat[12, 2] <- mean(Non_Changer$Morality) 
tp <- t.test(Changer$Morality,Non_Changer$Morality) 
tab6.mat[12, 3] <- ifelse(tp$p.value < 0.05, "Yes", "No")

tab6.mat[13, 1] <- mean(Changer$nuclear.knowledge) 
tab6.mat[13, 2] <- mean(Non_Changer$nuclear.knowledge) 
tp <- t.test(Changer$nuclear.knowledge,Non_Changer$nuclear.knowledge) 
tab6.mat[13, 3] <- ifelse(tp$p.value < 0.05, "Yes", "No")

tab6.mat[14, 1] <- mean(Changer$political.knowledge) 
tab6.mat[14, 2] <- mean(Non_Changer$political.knowledge) 
tp <- t.test(Changer$political.knowledge,Non_Changer$political.knowledge) 
tab6.mat[14, 3] <- ifelse(tp$p.value < 0.05, "Yes", "No")

tab6 <- data.frame("changer" = as.numeric(tab6.mat[,1]),
                   "non-changer" = as.numeric(tab6.mat[,2]),
                   "significance" = tab6.mat[,3])
xtable(tab6, digit=1)

###################################################
## Figure 2 
###################################################
fit5.1 <- glm(behavior.change ~ relevel(choice.at.treat, ref="No Change"), family = "binomial", data=df2019) ;summary(fit5.1)

p <- plot_model(fit5.1, type="pred", terms="choice.at.treat")+
  theme_fivethirtyeight()+
  theme(axis.text.x = element_text(size=10, face ="bold"))+
  theme(axis.title.x = element_text(size=10, vjust = -3, face ="bold"))+
  theme(axis.title = element_text(size=10, face ="bold"),
        title = element_text(size=12))+
  labs(title = "Predicteds Probabilities of Behavior Change by Attitude Strength")+
  theme(axis.title.y=element_blank())+
  theme(plot.title = element_text(hjust = 0.5, size = 20))+
  xlab("Attitude Strength")+ylab("Probability of Behavior Change")

pdf(file="Figure2.pdf", width=10, height=6, family="sans")
print(p)
dev.off()

###################################################
## Figure 3 
###################################################
#### 2019 Bar Graph 
Group <- df2019 %>% select(Treatment, opinion.change, behavior.change)
         names(Group)[2] <- "Opinion Change"
         names(Group)[3] <- "Behavior Change"

Group <- Group %>% pivot_longer(-Treatment, names_to = "variables", values_to = "value")
Group2019 <- summarySEwithin(Group, measurevar="value", withinvars=c("Treatment","variables"))
Group2019$variables <- factor(Group2019$variables, levels = c("Opinion Change", "Behavior Change"))

G1 <- ggplot(Group2019, aes(reorder(Treatment, -value), y=value, fill=as.factor(variables))) + 
  geom_bar(position=position_dodge(.9), colour="black", stat="identity")+
  geom_errorbar(position=position_dodge(.9), width=.25, aes(ymin=value-ci, ymax=value+ci)) +
  scale_fill_grey(start=0.9, end=0.2)+
  coord_cartesian(ylim = c(0,1))+
  scale_y_continuous(labels=scales::percent, breaks = c(0.1, 0.2,0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0))+
  theme_bw() +
  geom_hline(yintercept=0.5, linetype="dashed")+
  theme(strip.text.x = element_text(size = 13))

G2019 <- G1+labs(x="Treatment Group", y="% Opinion/Behavior Change", 
                 title = "Attitude Change Experiment (2019)")+
          theme_fivethirtyeight()+
          theme(plot.title = element_text(hjust = 0.5, size = 25))+
          theme(legend.position = c(0.50, 0.93))+
          theme(legend.title = element_text(size=15))+
          theme(legend.text = element_text(size=13))+
          theme(axis.text.x = element_text(size = 16))+
          scale_x_discrete(guide = guide_axis(n.dodge = 2))+
          theme(axis.title.x = element_text(size = 19, vjust = -1.5))+
          theme(axis.title.y = element_text(size = 15))+
          guides(fill=guide_legend(title="Type of Change"))

#### 2022 Bar Graph
## data1 <- df2022 %>% filter(Treatment!="Control") ## Data without Control group 

Group2 <- df2022.nocontrol %>% select(Treatment, behavior.change, opinion.change)
          names(Group2)[2] <- "Behavior Change"
          names(Group2)[3] <- "Opinion Change"

Group2 <- Group2 %>% pivot_longer(-Treatment, names_to = "variables", values_to = "value")
Group2022 <- summarySEwithin(Group2, measurevar="value", withinvars=c("Treatment","variables"))
Group2022$variables <- factor(Group2022$variables, levels = c("Opinion Change", "Behavior Change"))

G2 <- ggplot(Group2022, aes(x=reorder(Treatment, -value), y=value, fill=as.factor(variables)))+ 
  geom_bar(position=position_dodge(.9), colour="black", stat="identity")+
  geom_errorbar(position=position_dodge(.9), width=.25, aes(ymin=value-ci, ymax=value+ci)) +
  scale_fill_grey(start=0.9, end=0.2)+
  coord_cartesian(ylim = c(0,1))+
  scale_y_continuous(labels=scales::percent, breaks = c(0.1, 0.2,0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0))+
  theme_bw() +
  geom_hline(yintercept=0.5, linetype="dashed")

G2022 <- G2+labs(x="Treatment Group", y="% Opinion/Behavior Change", 
                 title = "Attitude Change Experiment (2022)")+
  theme_fivethirtyeight()+
  theme(plot.title = element_text(hjust = 0.5, size = 25))+
  theme(legend.position = c(0.50, 0.93))+
  theme(legend.title = element_text(size=15))+
  theme(legend.text = element_text(size=13))+
  theme(axis.text.x = element_text(size = 16))+
  scale_x_discrete(guide = guide_axis(n.dodge = 2))+
  theme(axis.title.x = element_text(size = 19, vjust = -1.5))+
  theme(axis.title.y = element_text(size = 15))+
  guides(fill=guide_legend(title="Type of Change"))

#### Merge Two Plots 

pdf(file="Figure3a.pdf", width=20, height=12, family="sans")
print(G2019)
dev.off()

pdf(file="Figure3b.pdf", width=20, height=12, family="sans")
print(G2022)
dev.off()


###################################################
## Figure 4 
###################################################
fit8 <- glm(opinion.change ~ nuclear.knowledge, data=df2019, family = "binomial") 
fit9 <- glm(behavior.change ~ nuclear.knowledge, data=df2019, family = "binomial") 

p1 <- plot_model(fit8, type="pred", terms="nuclear.knowledge")+
  theme_fivethirtyeight()+
  theme(axis.text.x = element_text(size=16, face ="bold"))+
  theme(axis.title.x = element_text(size=16, vjust = -3, face ="bold"))+
  theme(axis.title = element_text(size=16, face ="bold"),
        title = element_text(size=12))+
  labs(title = "Predicted Probabilities of Opinion Change \n by Nuclear Knowledge")+
  theme(axis.title.y=element_blank())+
  theme(plot.title = element_text(hjust = 0.5, size = 25))+
  xlab("Level of Nuclear Knowledge")+ylab("Probability of Opinion Change")

p2 <- plot_model(fit9, type="pred", terms="nuclear.knowledge")+
  theme_fivethirtyeight()+
  theme(axis.text.x = element_text(size=16, face ="bold"))+
  theme(axis.title.x = element_text(size=16, vjust = -3, face ="bold"))+
  theme(axis.title = element_text(size=16, face ="bold"),
        title = element_text(size=12))+
  labs(title = "Predicted Probabilities of Behavior Change \n by Nuclear Knowledge")+
  theme(axis.title.y=element_blank())+
  theme(plot.title = element_text(hjust = 0.5, size = 25))+
  xlab("Level of Nuclear Knowledge")+ylab("Probability of Behavior Change")

#### Merge Two Plots 
pdf(file="Figure4.pdf", width=20, height=10, family="sans")
grid.arrange(p1, p2, ncol=2)
dev.off()


end_time <- Sys.time()
print(end_time - start_time)
